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Abstract 

We introduce a mathematical model for ultrasound modulated optical tomog- 
raphy and present a simple reconstruction scheme for recovering the spatially 
varying optical absorption coefficient from scanning measurements with nar- 
rowly focused ultrasound signals. Computational results for this model show 
that the reconstruction of sharp features of the absorption coefficient is possi- 
ble. A formal linearization of the model leads to an equation with a Fredholm 
operator, which explains the stability observed in our numerical experiments. 
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1. Introduction 

During the last two decades, optical tomography (OT) has received signif- 
icant attention as a biomedical imaging modality. This can be attributed, in 
particular, to the fact that light at optical frequencies is harmless to the human 
organism and that optical properties of tissues reveal important biological infor- 
mation such as angiogenesis and hypermetabolism, both of which are well-known 
indicators of cancer [1] . Unfortunately, reconstruction in OT is also known to be 
severely ill-posed, and consequently the sharp imaging of optical properties is all 
but impossible. Various attempts to address this problem have been made. In 
this paper we are interested in a hybrid imaging method called Ultrasound Mod- 
ulated Optical Tomography (UOT, Q]) that combines the OT procedure with 
simultaneous modulation by a narrowly focused ultrasound beam in order to 
alleviate the instability of OT reconstructions. The idea is to combine the good 
tumor specificity of OT with the high spatial resolution of ultrasound imaging. 
This approach utilizes the experimentally observed interaction between ultra- 
sound and light propagation in tissue [H [J. In UOT, a coherent light source 
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irradiates the tissue sample and causes interference patterns to form on the 
surface of the object, so-called speckles. A narrowly focused ultrasound wave 
is simultaneously induced in the tissue, influencing its optical properties and 
thus modulating the speckle pattern with ultrasound frequency. By measuring 
properties of this modulation, information about the incident light intensity at 
the focus location of the ultrasound beam can be obtained. Hence, by scanning 
the focus of the ultrasound wave throughout the sample, a quantity related to 
the light intensity in the object's interior can be determined. This type of inter- 
nal information is usually not available from OT measurements due to multiple 
scattering of photons in optically dense media, although there are other vari- 
ants of optical tomography that also strive to recover this information (e.g. [3]). 
It can be expected that this additional knowledge can help in stabilizing the 
inversion process and render it substantially less ill-posed than the original OT 
problem. For the UOT model we present in this paper, numerical experiments 
and an initial analysis suggest that this intuition is justified. 

The literature contains a number of models that address the UOT tech- 
nique, see for example [H HI [51 03 H] • Most of them describe the coupling 
between ultrasound and light in terms of stochastic quantities, which permits 
particle-based simulations of the light intensity modulation effect caused by the 
ultrasound wave. On the other hand, for optical imaging in turbid media at 
a depth of several centimeters, photon intensities can be accurately modeled 
by the diffusion limit. Under certain assumptions, this allows us to formulate 
a model for the UOT procedure based on a parameter identification problem 
for a set of coupled diffusion-type partial differential equations. This model, 
along with a description of the measurements is presented in Section [2j In 
Section [3j we outline a simple algorithm that can be used to reconstruct the 
spatially varying absorption coefficient from UOT measurements with focused 
ultrasound signals. Examples of the resulting reconstructions for numerical 
phantoms are provided in Section [4j In Section [5j we formally linearize our 
model and obtain an equation that relates perturbations in the absorption co- 
efficient to those in the measurements by a Fredholm operator acting between 
appropriate Sobolev spaces. This provides a partial explanation to the stable 
reconstruction observed in our numerical experiments. The last section contains 
final remarks and conclusions. 

2. Mathematical model 

A detailed description of the physical underpinnings of the UOT procedure 
can be found, for instance, in [TJ Ch. 13]. We give a brief description of the 
set-up here. 

Let the object of interest occupy the domain C 1R 3 . The internal optical 
properties in the diffusion limit are described by the reduced scattering coeffi- 
cient fj,' s and the absorption coefficient fi a . For imaging soft tissues, it is common 
to assume fj,' s roughly equal to a known constant throughout Q, while the spa- 
tially varying absorption fi a {x), x 6 f2, represents the target of reconstruction. 
It is also assumed that the tissue of interest is turbid (highly scattering) , so that 
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Ha{x) <C fi' s . It is known that in such media, the light intensity u(x) inside fl 
can be accurately described by the diffusion approximation (e.g., [HE]). 

It has been shown experimentally that coherent light can be modulated by 
an ultrasound field inside the turbid medium 2J. Various explanations have 
been put forward for this effect [5]. 

The experimental setup in UOT involves dealing with the time dependent 
light intensity of individual speckles. The model presented below is derived 
under two assumptions, which are satisfied in standard UOT applications [1]: 

• Weak scattering assumption: The optical wavelength is much shorter than 
the mean free path. 

• Weak ultrasound modulation assumption: The ultrasound-induced change 
in the optical path length is much less than the optical wavelength. 

The measured signal is the autocorrelation function [2] at a detector location 

Gifa.r) = (E(r],t + T)E*(r),t)) t , 

where angle brackets denote averaging over time, and the electric field E is 
related to the light intensity I as I(i],t) — \E(r) 7 t)\ 2 . It has been shown ex- 
perimentally [5] that over time scales r » lfis coherence of the exiting light is 
lost, i.e. Gi(ji,t) — > as r — > oo, due to the Brownian motion of scatterers. 
However, on short time scales on the order of the period of the ultrasound field 
- i.e. the regime we are interested in -, Gi(t],t) has been observed to oscillate 
at the ultrasound frequency. We will therefore neglect contributions from the 
Brownian motion of scatterers since it is unrelated to the ultrasound field. In 
the absence of an ultrasound field, and on these time scales, we would then have 
Gi(rj, r) = const. In the following, we will derive expressions for G\ and, in 
particular, its modulation depth, i.e. the magnitude of the oscillation of G\ at 
the ultrasound frequency. We will then relate these quantities to solutions of 
partial differential equations that we will use for our reconstruction scheme. 

A path integral model. For a point source of unit strength at a location a, and 
a detector measuring photons exiting the domain at r\ S dfl, we can write 

G 1 (a,r 1 ,r)=P d G(a,r ] ,T), G(a, V> r) = £ P s (E s (t + r)E* s (t)) t 

s—s(a 1 ri) 

where the sum extends over all paths s that connect source a and boundary 
location rj. P s is the fraction of the incident intensity that scatters along s 
multiplied by the probability of a photon not getting absorbed along this path. 
P d is the probability that a photon that makes it to a point r\ on the boundary 
is able to cross the boundary from tissue into the detector. E s then denotes 
the phase of the electric field at rj of photons following path s. Consequently, 

(E s (t)E* s (t))t = 1- 

Consider now the situation in which the ultrasound field p(x, t) induces phase 
shifts d(j)(x,t) on all paths along an infinitesimal path element ds(x). As shown 
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in [5] , such phase shifts can be induced both by the periodic motion of scatterers 
in the ultrasound field as well as by the modulation of the index of refraction 
by the pressure field. We then have 

(E s (t)E* s (t + T)) t = (e-p(-ij[^^) 

where integrals are assumed to be along a path s from a to rj. By computing 
how the index of refraction and the phase shifts induced by scatterer movement 
depend on an ultrasound pressure field with frequency cu a , we can use the results 
in [5] to write above expression as 

(E s (t)E* s (t + T)) t = ^f s exp[-a|p(z)| 2 (l- coster)] ds, 

where a is a proportionality constant and \s\ is the length of path s. (Note in 
particular that the proportionality to the square of the pressure has also been 
observed experimentally, see [2].) Consequently, 

Gi(<t,77,t) = P 9 ^A^exp[-a|p(ir)| 2 (l-cosw a T)] ds. 

S=s(<7,7]) 

As has been shown experimentally [2], the temporal variation of the exponent 
is relatively small. We can therefore approximate 

l-^j>(*)| 2 (l-COS WQ T)d S ] . (1) 

It follows that we can write the autocorrelation function as the sum of two 
terms: 

G 1 (a, V , t) = G 1 (a, r,, 0) - aP 9 E S = S ( CT)J7 ) P. ft J. \p(x)\ 2 (l ~ cosc Q r) ds. 

The first of these is the time average light intensity, whereas the second is the 
temporal variation of the autocorrelation function due to the ultrasound field. 
To first order in the small parameter a, this expression equals 

G x (cr, rj, t) = Gi {a, r,, 0) - aP 9 f n G(a, x, 0) \p(x) \ 2 G{x, r), 0) dx (1 - cos uj a r) . 

Finally, if light is incident with an intensity S(a) at source positions a G dfl, 
the overall autocorrelation function at detector location 77 can be written as 

Gi(v,r) = / an 5(<7)Gi(a,77,r) da. (2) 

Using the previous equation, and defining the time averaged light intensity 
u(x) = J an S(<j)G(cr, x, 0) da for all x G O U dfl, we can then write 

Gi(t),t) = P d u{rj) - aP 9 [ u(x)\p{x)\ 2 G{x,ri,0) dx (1 - cosw a r) (3) 

Jn 

= P a [ M (77)- W (77)(l-cosw a r)], (4) 



G^a^r) =P 9 E s = s{ „, v) Ps 
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where 

v(r/) = a / u(x)\p(x)\ 2 G(x, rj, 0) dx. (5) 
Jq 

This representation of the correlation as a sum of a time averaged photon flux 
plus a temporally variable term has given rise to the name tagged photons to 
denote v(x). Equation ^ makes it clear that tagged photons originate at the 
site x of interaction of the steady-state light field u(x) and ultrasound field p(x). 
However, since Gi(?7,t) is not a photon flux but a correlation function, we will 
not use this term any further. 

In our reconstruction algorithm below, we will assume that the amplitude 
v{rf) of the temporal variation of Gi(?7,t) - i.e. the modulation depth - is the 
measured signal. While the time average u(rj) can also be measured, using it 
for inversion leads to the diffuse optical tomography problem that is known to 
be severely ill-posed. 

A partial differential equation model. For our reconstruction algorithms, we 
would like to relate our signal v(rj) to the solution of a partial differential equa- 
tion. To this end, note that G(x, y, 0) = J2 s=s (x y ) ^ * s tne time avera g e prob- 
ability that a photon starting at x is found at y. For the turbid medium that 
we consider in this contribution, light propagation can be accurately described 
by the diffusion approximation in which photons perform a random walk. The 
time averaged light intensity u(rf) = S(a)G(a, 77, 0) da must then satisfy the 
following equation: 

-V -DVu(x) + n a (x)u(x) =0 in fl, (6) 

where 

is the diffusion coefficient. Due to the assumptions stated at the beginning of this 
section, D w ^jp- w const. To simplify the notation, we set /1 := fi a in the rest 

of the text. Equation ^ needs to be completed by boundary conditions. For 
tissue in contact with a surrounding medium, Robin-type boundary conditions 
are typically chosen [9]: 

2D 9 ^ + War) = S(x) on dQ. (8) 

on 

Here n denotes the outward normal to the surface 9f2 and 7 > is a constant 
describing the optical refractive index mismatch at the boundary, and is related 
to P s . In particular, the assumptions underlying the diffusion approximation 
imply that P d < \ and < 7 < 1, and 7 = 1 if P d = \. 

On the other hand, to represent v(r\) as the solution of a partial differential 
equation, we have to consider the equation that G satisfies. G(x,y,0) is the 
probability that a photon originating at x reaches y, absent an ultrasound field. 
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For random walk models, it is known that G(x, y, 0) satisfies a diffusion equation 
[TUl ITT] , which in our case is 

-V-DVG(x,y,0) + fi(x)G(x,y,0) = S(x-y) in Q. 

The question of boundary conditions is less clear. It is well known that if every 
particle that reaches the boundary leaves the domain (i.e. P d — 1), then the 
correct boundary condition to choose is G\gn — 0. On the other hand, if all 
photons are reflected and none can leave (i.e. P a — 0), then n- VG|asi = is the 
correct boundary condition. In either of these two cases, n ■ VG|an is the flux of 
particles across the boundary. However, we have been unable to find literature 
on the case < P d < k (see, however, [12] for the case where each particle that 
reaches the boundary is replaced by more than one new particle, a situation that 
formally corresponds to the situation where the fraction of particles that can 
leave the domain satisfies P d < 0). Since intuitively, G denotes a photon flux, 
we conjecture by way of analogy that G also satisfies Robin boundary conditions 

2D dG(x y,0) + Q, y ) =0 on d Q_ ( 9 ) 
on 

Under this assumption, we have that the amplitude v(rj) (up to the constant 
factor P d ) of the time variation of the autocorrelation function Gi(r/, r) satisfies 
the following boundary value problem: 



—V • DVu(i) + n(x)v{x) — a\p(x)\ 2 u(x) in ft, 
2D^+jv(x)=0 on dft 



Note that if we were to view v as a fluence of virtual or tagged photons, then 
our conjecture implies that the equation for this virtual fluence has the same 
boundary conditions as that for the incident fluence u. 

Measurements. In principle, the interferometric detectors for the modulation 
P d v(rj) visible beyond the boundary could be placed along the entire boundary. 
In practice, however, we will only be able to measure at a small number of 
locations. To simplify the discussion, we will assume in the following that only 
a single detector is used. More elaborate experimental setups could use multiple 
detectors to suppress the effects of noise on the reconstruction. 

The inverse problem. We can now formulate the inverse problem addressed in 
this work: Assuming that for a given point n G dQ and a number of ultrasound 
fields p^(x) indexed by £, the values 

h(0--=vHv) (ii) 

are known in the coupled system of equations 

-V • DS7u(x) + fi(x)u(x) = in Q, 

2D d ^+ 7 u(x) = S(x) on dn, 

-V • DVvt(x) +(jt(x)v*(x) = a\pt{x)\ 2 u{x) in ft, [ > 

2Z>^M + 7 u«(a;) = on dn. 
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Then recover the absorption coefficient fi inside a region of interest U C £1 with 

u c n. 

We remark that in the applications of ultrasound modulated optical tomog- 
raphy available in the literature, the ultrasound pressure field p{x) is always a 
beam focused on a single point. In particular, the algorithm we show below 
is based on the assumption of perfectly focused beams |p^(a;)| 2 = S(x — £), al- 
though we will test in Section [6] how the algorithm performs on data for which 
this assumption is not satisfied. The formulation above is more general in that 
it allows arbitrary fields p(x). An application of this includes ultrasound pres- 
sure fields that are focused not on points but on spherical surfaces for synthetic 
focusing, as mentioned in Section [6j 



3. Reconstruction algorithm 

In this section, we introduce a simple algorithm that can be used to compute 
numerical reconstructions for the above inverse problem. In the following, we 
will make the assumption that the pressure field is perfectly focused on a location 
£ G Q, i.e. |p£(x)| 2 = 5{x — £). As discussed in Section [6j this is of course not 
practically feasible, so our assumption is understood to mean that the real 
pressure field approximates a perfectly focused one. 

Let G(x,y) be the Green's function for the diffusion model i.e. the 
solution of 

( -V x -DV x G(x,y)+n(x)G(x,y)=5(x-y) x G ft, 

\ 2D^l +1 G(x, y)=0 xedQ. [l6> 



Then, ( 12 ) implies 

vHx) = aG(x,Ou(0, 

and thus, 

h(0 = aG( V ,Ou(0, u(0 = 



MO 
aG(n,0' 



Substituting this expression for u into the first equation of (12), we obtain an 
equation for recovering fi: 

The apparent difficulty in using this formula for reconstruction is that it is 
implicit in \l since both D and the Green's function G depend on the absorption. 



However, we can construct the following natural iterative scheme for ( 14 ) : 



Initial step: Using an initial guess /jP for the absorption coefficient (e.g. 
fi = const), compute the corresponding Green's function numerically, and 



apply formula ( 14 1 to find a new approximation for the absorption. 



Iterative step: Using the current approximation /i fc , re-compute Green's 



function and D and apply formula ( 14 ) to find an updated absorption 
coefficient fi +1 . 
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We do not consider the convergence properties of this scheme here, but note that 
in our numerical tests presented below the iterates converged reliably, albeit not 
very rapidly. 



Numerical implementation 

Implementation of the algorithm outlined above requires the following steps: 
• Simulation of the forward model to generate synthetic measurements, 



repeated computation of the Green's function G(x,y) for equation (13 1, 



repeated evaluation of the iteration formula ( 14 1 



These steps are discussed in the following subsections. In this work, we only con- 
sider measurements obtained by forward calculations from mathematical phan- 
toms, rather than actual experimental data. All computations were done in 2D, 
although they can be readily carried over to 3D. For the finite element calcu- 
lations involved in the reconstruction scheme, the Open Source finite element 
library deal. II [13l [14] was used. 



4--1- Forward simulations 

In order to generate the measurements h(£) (s ee |TT| ), we need to compute 
the solution u(x),v^(x) of the forward problem (|12[) for a set of given data 
D, /i, S (diffusion coefficient, absorption coefficient, incoming light flux) and 
an ultrasound signal focused at the point £ 6 U. Then, evaluating at the 
detector location ry, we obtain the measurement value h(£). 

4- 1-1. Computational setting. 

We take ft to be the square [0,5cm] 2 , which approximately corresponds to 
the relevant dimensions in practical applications. For the boundary light source 
S in @, dQ is split into dflx = {x £ dfl : x x = 0} and dtt 2 = dfl\ 0fi x . 
Constant illumination is assumed on dfli and no photons are injected on dD,2'- 

The modulation depth is measured at a single detector location n — (5cm, 2.5cm). 
This layout is depicted in Fig. [T] 

4-. 1.2. Incident light field. 

Since in our model the incident light intensity u is independent of the shape 
and location of the ultrasound waves in the tissue, u only needs to be computed 
once. For this computation, a finite element approximation to u is constructed 
on a regular rectangular grid using Qi finite elements |15j . solving equations 
|6| ([8]) . The left panel of Fig. [2] shows u for the case of a constant absorption 
coefficient fi. 
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Figure 1: Setting of numerical experiments: Domain O, area of interest U, incident light 
source S(x) on the left, and detector point r\ on the right. 



4-1-3. Ultrasound field. 

In our numerical examples, we use Gaussian-shaped synthetic ultrasound 
signals: 

p(x) = Cexp(-J2 l -^Pj, (16) 

where C is a normalization constant. By choosing different variances ex 2 we can 
model varying focusing properties of such pressure field. 

To simulate scanning of the ultrasound focus, focusing points i — 

1, . . . , N} are placed at the vertices of a square grid covering the area of in- 
terest, here chosen as the square U = [0.5cm, 4.5cm] 2 C fi. For each i we then 
construct a signal p^ (x) focused at fj by setting 

p^(x) :=p(x - &). 

To simplify notation we set v l := and p l := pfi . 



4-. 1-4- Modulated light field and measurements. 

Given u and \p l \ 2 , we compute the intensity of the modulated light v l (x), 
using equations |l2]). The equations are again solved using Qi finite elements. 
Two examples for v l are shown in Fig. [2] for two different focus positions. The 
modulated light intensities v l arc then evaluated at the sensor location r\ to 
yield the measurements h(£i) — v l (rj). 



4-2. Green's function and reconstruction 

The reconstruction algorithm requires knowledge of the Green's function G, 
which, given the absorption coefficient fi and resulting diffusion coefficient D, 
solves (13). Hence, we compute G by solving another diffusion problem with 
homogeneous Robin boundary conditions and a suitable approximation to the 
delta function on the right hand side. As before, this is done using a finite 
element scheme, where we choose a different, coarser mesh than in forward 
problem calculations to avoid committing inverse crimes. 
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Figure 2: Left: Incident light intensity u for constant absorption coefficient. Center and right: 
Modulated light intensity for two different focus points £. Note that v depends on the focus 
position as well as the intensity of u at the focus. 



An obvious problem in the reconstruction formula ( 14 ) is that it involves 
derivatives of the measurement data h(£), which causes instabilities in the pres- 
ence of noise. Possible regularizations for this problem are well-studied (e.g. 
[5]), and the stability analysis in Section [5] suggests that this is the only source 
of instability in the reconstruction process. Hence, we opt not to add extra reg- 
ularization and compute the derivatives by a simple central finite differencing 
scheme. Without adding noise to the measurements, it turned out that in all 
of our computational experiments, the regularization stemming from discretiza- 
tion on a fixed grid was sufficient for convergence of the iterative scheme based 



on (14). 



If.,3. Numerical phantoms 

To test our algorithms, we use three test cases in which the true absorption 
coefficients have the following form: 



A disk-shaped inclusion K C fi with midpoint (2.5cm, 2.5cm) and radius 
0.5cm. The absorption coefficient is assumed to be equal to /2 outside the 
inclusion and slightly higher inside: 



(i , x e fi \ K 

1.2/2, xEK. 



For the same inclusion K, a much higher absorption coefficient contrast 
lf{x) -- 



[i , x e \ K 
10 p, , x e K . 

A more complicated coefficient with multiple inclusions of different mag- 
nitude between 1.2/2 and 2.0/2. Their exact shape is shown in Fig. [3] This 
case tests the ability of our algorithms to resolve several nearby objects. 

For actual numerical values, we used /2 = 0.023cm -1 , /i' s = 10.74cm" 1 and 
7 = 0.431cm" 1 in our computations. These values represent typical optical 
properties of soft tissue [16] . 
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32 - — -0.500 — 0.070 

29 ^^^^^"^^ 

^-0.00 ^BB^B^ *0.0 

Figure 3: Test cases for absorption coefficient fi* . 

4-. 4- Reconstruction results 

For the results shown in this section, measurements were produced using 
the ultrasound signal arising from setting variances a\ — o~2 — 0.1cm in the 
Gaussian ( |16[ ), resulting in sharp focusing in each direction (see the center panel 
of Fig. [5] below) . Fig. [4] shows reconstructions of the three different absorption 
coefficients for scanning the ultrasound focus on a 100 x 100 mesh of points 
inside the area of interest U. 

p-0.032 ,- 0.500 ,-0.070 ^ 

-0.02, ^^^^^ 0.375 ^^j^^ ^A^^ 

0.026 ^/ BJP^ 1 0.250 J y*"* I 

^BBJJJJJHB^r ^BBJ BJr 

^-0.020 ^^^^^^^^ ^-0.0 

m 0.032 H 0.500 m 0.070 

N = 40 I N = 70 I jf\ N = 41 

-0.029 / /'~'~* ^^k. .^BP^^^H^ 0.053 l^^l 

-0.026 I 0.250 ^fl ^ I 0.035 ^^^^H 

^BBJ| ^BBJ B-0.018 ^BBJ / 

^-0.020 y ^-0.00 ^-0.0 



Figure 4: Reconstruction results for the three coefficient cases: after the first step of the 
algorithm (top) and after N = 40, 70 and 40 iterations, respectively (bottom). 

The principal observation from these results is that under the main assump- 
tions of the model, i.e. turbid medium (and thus fi <C (i' s ), virtual light source, 
and strong focusing, our reconstruction scheme has four desirable properties: 

• It converges, even for the second case where (i) we start far away from the 
exact coefficient and (ii) the exact coefficient has a large dynamic range. 

• It is stable, i.e. the errors introduced through discretization of the equa- 
tions, finite differencing of data, and using different meshes for reconstruc- 



40 



11 



tion and generation of synthetic data do not lead to inaccurate reconstruc- 
tions. 



• It can recover sharp interfaces without excessive blurring. 

• It can recover quantitatively correct values of absorption. 

These are significant advantages compared to many other optical tomographic 
methodologies. 

5. Stability of the linearized problem 

The quality of reconstructions shown above, especially the recovery of sharp 
singularities, is at first surprising, given that the standard OT problem is 
strongly ill-posed. In this section, we will make a first step towards under- 
standing the stability of the UOT procedure. 



Note that even though equations ( 12 ) defining u and v are linear, the relation 



between the absorption coefficient /i and the measurements h is nonlinear. In 



this section, we consider a (formal) linearization of the system ( 12 1 that will 
allow us to gain some insight into the local properties of the inverse problem. 

Let n C R d with d = 2 or d = 3 be an open bounded domain with C - 
boundary. We use a formal linearization, assuming that /i is a small perturbation 
of a known absorption /io > 0,/io € C 0,1 (f2), and then applying the formal 
asymptotic expansions 

H(x) = no(x) + £fii(x) + o(e), 
u(x) = u (x) + eui(x) + o(e), 
v^(x) — Vq(x) + evf (x) + o(e), 

where e — > 0. Our goal is to relate the first order perturbations of the absorption 
coefficient fx% and the measurements := vf(rj), where rj £ 917 is the location 

of the detector. 

Let us again assume perfectly focused ultrasound, i.e. |p^(x)| 2 = S(x— £). By 



inserting the above expansions into equations ( 12 ) and sorting terms according 



to powers of e, we then get the zeroth order perturbation system 

- V ■ DVu Q (x) + t i a (x)u Q {x) = 0, (17) 

-V ■ DVv%(x) + fio(x)v§(x) = a6(x- £)ua(x), (18) 

and the first order perturbation system 

- V • DVui(x) + /j,o(x)ux(x) — —fj,i(x)uo(x), (19) 

-V • DVv\(x) + fi (x)vl(x) = a5(x - g)ui(x) - m(x)v§(x) (20) 

for all x G CI, complemented by inhomogeneous Robin boundary conditions as in 

([8| for mq and homogeneous Robin boundary conditions for Vq, u\ and v\. Here 
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we neglected the (weak) dependence of D on fj, and instead set D = const > 
for the rest of this section. 

Equations (fl7|)-( 18 ) imply that uo and v$ are solutions to the forward model 
for absorption coefficient fiQ. The standard elliptic regularity theorems (e.g., 
[T7| ) imply u G H 3 (il), and by the Sobolev embedding theorem u G C 1 (J7) 

Let us assume that the absorption coefficient is known near the boundary, 
so that it suffices to consider perturbations fi\ supported in an open set U with 
C 2 -boundary such that U C f2. We assume the data to be given for all 

£ G {/. In what follows, we derive an explicit formula for the dependence of /ii 
on /ii and then study properties of the corresponding linear operator. 

Let us denote by Go(x, y) the Green's function as defined in ( 13 ) correspond- 
ing to the background absorption coefficient /iq. Equation (18) implies that for 
all x G fl and (£[/, 



aGo(x, z)5(z — £)uq(z) dz 
uG {x,£,)u (O- 



From ( 20 ) we can now deduce that 



aG (a;,£)ui(£) - « u o(C) / G (x, z)G (z, (z) dz. 



Evaluating at x — r\ and solving for u\ yields 
hi(0 uo(0 



aG (?7,C) G (r),£) J n 



G (v,z)G (z,()ix 1 (z) dz. 



We now use this expression to eliminate u\ from (19). Noting that the differen- 
tial operators now act on £ and that 



DV £ + fx (0}G (x,0 = S(x-0, 
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we get 



= «o(OA*i(0 + [-Vf-DV € +Ai (0] 



-V 6 -DV 6 + MO] 



«o(0 



G (77,0 y n 
«o(0^i(0 + [-Vr-DV { + /*o(0] 

«o(0 



Go(r7,z)G„( "-') 

MO 







- 2L> 




"0 


(0 



Go (77,6 



Go(»?,0 
VP 
Go(t?,0 

G (r?,e)Mi(O 



G (?7, z)G (z,£)/ii(z) dz 
/ G (»7, z)G (z,C)a i i( z ) ^ 



We will frequently view Gq(j], y) as a function of y in the following and hence 
introduce the notation 

G v (y) :=Go(r7,y) for y G [T. 

Note that since 77 € 90, Gq has no singularities on U and hence is a regular 
solution to (17 1 there. The elliptic regularity and Sobolev embeddings imply 
Gg G C x (jJ). 

Let us define the following operators acting on functions g defined on U : 
1 



Kig(0 := - 
K 2 g(0 
and 



2u {0 
D 

MO 



[-v« • DVd 

MO 
Gg(0 



\M0~ 




[g v o(0. 


u 



G%(z)Go(z,09(z) dz, (21) 



Go(z)G (z,Og(z) dz 



(22) 



F :=1 — K\ — K 2 . 



In terms of these operators, our considerations above imply that /zi is a solution 
to the following linear equation: 



FMO 



1 



2uq(0 



-v r Dv e + M0] 



aG%(0 



(23) 



In order for the above expressions to be well-defined, we have to make sure 
that ito and Gq are bounded away from zero on U. The following lemma follows 
immediately from the Hopf Lemma (e.g., [I9j [20]): 

Lemma 1. There is a constant c > such that uo > c and Gq > c on U . 

Next we consider the properties of the integral term involved in K± and K%. 
The important observation here is the following: 
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Lemma 2. The mapping 



H- / G (z,-)G^(z)g(z)dz (24) 



is a bounded linear operator from L 2 (U) to H 2 (U). 

Proof: Let us assume that g £ L 2 (U). Since Gq € C(U), multiplication by 
Gq is a bounded linear operator on L 2 (U). The following integration against 
Gq(z, •) results in the solution to the diffusion equation with homogeneous Robin 
boundary condition and right hand side G^g £ L 2 (U). Elliptic regularity theory 
(e.g., [13 [H]) implies that this is a continuous operator from L 2 (U) into H 2 (U). 
□ 

Because of the compact embedding of H 2 (U) in L 2 (U), the operator defined 



by (24), viewed as a mapping from L (U) to L (U), is compact. In (21), this 

MO 



operator is multiplied by the factor 

1 



2«o(0 



[-v« ■ Z?V £ ] 



Gg(0 



(25) 



The functions u ,Vwo,G[{ and VGq are all bounded on U because uq,Gq S 



C 1 ([/). Since m and Gq satisfy (17), the terms V • D\7u and V ■ -DVGq are 
bounded on U as well, and u^ 1 and (Gq) -1 are bounded due to Lemma [lj 



Consequently, multiplication by (25) represents a bounded linear operation on 
L 2 (U), and so K\ is a compact operator in L 2 (U). Similarly, K2 is a compact 
operator in L 2 (U). This leads us to the main result of this section: 

Theorem 3. F : L 2 (U) — > L 2 (U) is a Fredholm operator of index zero. 

Thus, the kernel Af(F) of F has finite dimension and the range 71(F) is closed 
and of finite codimension, equal to the dimension of the kernel. This immedi- 
ately implies the following result: 

Corollary 4. F as an operator from the quotient space L 2 (U)/Af(F) to 1Z(F) 
has bounded inverse, and the following norm equivalence holds: 

cill-F/Hz^cto < ll/IU 2 (c/)/AA(F) < c 2 \\Ff\\ L 2( U) . (26) 



The L 2 -norm of the right hand side expression in ( 23 1 can be estimated in terms 
of the ff 2 -norm of the measured perturbation hi , so that we obtain the following 
stability result: 

Theorem 5. Under the stated assumptions, there is a constant C > such 
that the following relation holds: 

\\ih\\li(u)/n-(f) < CWhiWxnuy (27) 
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We conjecture that the kernel Af(F) is in fact trivial, and thus the operator F 
is invertible. This would imply that fii is uniquely determined by the measured 
perturbation hi, and allow us to replace the quotient space norms in (26) and 
(27) with the regular L 2 norms. However, we have not been able to prove this 
result yet. 

Smoother norm coercive estimates for the absorption can be obtained if more 
is assumed about the unperturbed absorption /j,q and the domain. For instance, 
if no € C°°(fi), S G C°°(dil), and fi has smooth boundary, the operators Ki 



and K 2 defined in (|21j)-([22j), are of order —2 and —1, respectively, in the Sobolev 
scale: 

Ki : H S {U) -> H S+2 {U), 
K 2 : H"(U) -> H S+1 {U). 

This and the Sobolev embedding theorem [3T] imply that for any s > 0, F is 
Fredholm as an operator 

F : H S (U) -> H S {U). 
This, in turn, leads to the estimate 

< c(\\f\\ LHu) + \\Ff\\ H s {u) ) (28) 
for all / £ H S (U). Thus, we have the following result: 

Theorem 6. Under the stated assumptions, for any s > there is a constant 
C such that 

\\Hi\\hs(U) < C (llMl|U 2 (t/) + II^i||h=+ 2 ((7)) • 

Clearly, if only a specific value of s is of interest, the smoothness assumptions 
on /io, S and dVt can be relaxed appropriately. 



6. Conclusion and outlook 

In this paper, we have introduced a partial differential equation model of 
ultrasound modulated optical tomography to derive a simple reconstruction 
scheme for recovering the spatially varying absorption coefficient from boundary 
measurements. While we could demonstrate stable, sharp and quantitatively 
accurate reconstructions, some of the assumptions made here need to or can be 
improved upon for practical applications. In particular, these are: 

Detector locations. In the discussion of stability above, as well as in our numer- 
ical reconstructions, we have chosen a single detector point 77. However, using 
detectors distributed over a part T of the boundary dfl should help to suppress 
the effect of noise in the measured data. 
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Ultrasound signal with elongated focus. In practice, perfect focusing of ultra- 
sound waves is not a realistic assumption [22 . How well an ultrasound signal 
can be focused depends, in particular, on the geometry and bandwidth of the 
transducer. For example, it is known from experimental measurements (e.g., 
[53]) that focused ultrasound signals have an intensity profile similar to the one 
shown in Fig. [5] (left). This signal has significantly sharper focus in the direction 
transverse to the transducer lens, while the well-focused Gaussian signal used 
in our results does not reflect this behavior. 



Figure 5: Left: Simulated ultrasound pressure field \p\ 2 with transducer at the bottom. Mid- 
dle: Gaussian ultrasound signal \p\ 2 with <ri = <72 = 0.1. Right: Gaussian signal with 
(71 = 0.1, CT2 = 0.3. 

To illustrate the effect of relaxing the assumption of perfect focus, we com- 
puted reconstructions for the case where the ultrasound intensity is a Gaussian 
signal with sharp focus in x-direction and elongated focus in y-direction (Fig. [5j 
right). As in the previous section, the ultrasound focus is scanned on a 
100 x 100 mesh to produce synthetic measurements. At the same time, the 
reconstruction algorithm is left unchanged, i.e. still assumes perfect focus. 

Reconstruction results are shown in Fig. [6] The deterioration of the recon- 
struction - in particular in the direction of the ultrasound beam - is obvious. 
The results also contain artifacts at the vertical boundaries and close to the 
detector location. A more sophisticated reconstruction scheme might be needed 
to treat the non-perfect focusing in these calculations. 

Synthetic focusing. Instead of attempting to perfectly focus the ultrasound 
waves in space, synthetic focusing allows the use of non-localized ultrasound 
fields and reconstructs the signal by superposition. This approach was sug- 
gested in [23]: It combines various basis sets of non- focused ultrasound waves 
(e.g., spherical or monochromatic planar ones), with a post-processing step that 
synthesizes the would-be response to a focused illumination. In particular, in 
the case of spherical waves, the post-processing (synthetic focusing procedure) 
is essentially equivalent to thermoacoustic tomography inversion (see |25j). We 
plan to investigate the applicability of this approach to UOT in the future. 

Uniqueness of reconstruction. Proving uniqueness of reconstruction, both in 
the non-linear and linearized versions, still remains a challenge. In particular, 



we conjecture that the operator F in (23 1 is in fact invertible, and thus there 



is uniqueness of solution of the linearized problem, which would replace the 
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Figure 6: Reconstruction results for ultrasound signal with elongated focus: after the first 
step of the algorithm (top) and after N iterations (bottom). 



quotient space norms in (26 1 and (27 1 with the regular L norms. At the same 



time, a complete characterization of the kernel of the operator F is non-trivial 
and left for future work. 

Summary. Despite these opportunities for future work, in this paper, a diffusion 
based model is provided for the ultrasound modulated optical tomography pro- 
cedure using well focused ultrasound waves. An iterative algorithm is suggested 
to recover absorption from measurements of the amplitude of ultrasound mod- 
ulation. The provided numerical results show feasibility of the algorithm and 
possibility of good reconstructions, both with regard to locating sharp interfaces, 
as well as recovering correct numerical values of the absorption coefficient. Such 
stability and resolution are impossible to achieve in standard optical tomogra- 
phy. The stability of reconstructions is explained by the stability estimates 
derived in Theorems [5] and [6] for a linearized model. 
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